
## PEASANT RESISTANCE AND ECONOMIC AFFLUENCE: LESSONS FROM PARAGUAY
## Appendix B: Robustness Tests
## January 2024

## Code Authors:
## Germán Feierherd (gfeierherd@udesa.edu.ar) 
## Jorge Mangonnet (jorge.mangonnet@vanderbilt.edu)


#### Preamble ------------------------------------------------------------------

## Workspace setup
rm(list = ls())		         # clear all objects in memory
options(max.print = 5000)  # expand display
options(scipen = 999)      # disable sci notation
dev.off()                  # reload graphic device
cat("\014")                # clear console

## Load/install packages
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  readxl,
  interflex,
  marginaleffects,
  haven,
  readr,
  reshape2,
  magrittr,
  dplyr,
  plm,
  lfe,
  pcse,
  lmtest,
  stargazer,
  xtable,
  texreg,
  estimatr, 
  fixest,
  margins, 
  sjPlot, 
  multiwayvcov,
  #  rgdal,
  ggplot2,
  wesanderson, 
  fixest,
  did,
  sf, 
  sandwich,
  AER,
  MASS)

## Set working directory
gf <- "~/Dropbox/Working_Papers/investigacion protesta Paraguay/"; main <- gf  # german's path
jm <- "~/Dropbox/Research/investigacion protesta Paraguay/"; main <- jm        # jorge's path
setwd(main)

## Functions
source(paste0(main, "4_code/functions_paraguay.R"))

## Load full dataset
load(paste0(main, "4_code/DataSetPY.Rda"))


#### Appendix B: Robustness Tests ---------------------------------------------

#### Table B.1: Assessing Reporting Biases ----

# Munis w/ zero resistance events
zeros <- dta %>%
  group_by(codigo) %>%
  summarise(n_conflict = sum(n_conflict)) %>%
  dplyr::filter(n_conflict == 0) %>%
  dplyr::select(codigo)
# Define restricted sample
dta_r <- dta %>%
  subset(!(codigo %in% zeros$codigo))

# Model 1: Probit, w/o controls
m1 <- glm(b_conflict ~ I(ln_suitability*ln_price) + factor(codigo) + factor(year), binomial(link = "probit"), data = dta); summary(m1)
se1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC2", cluster = "group")); se1

# Model 2: Probit, w/ controls
m2 <- glm(b_conflict ~ I(ln_suitability*ln_price) + ln_pop + factor(codigo) + factor(year)*distance_hubs + factor(year)*untitled_share, binomial(link = "probit"), data = dta)
se2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC2", cluster = "group")); se2

# Model 3: Restricted Sample, w/o controls
m3 <- plm(ln_conflict ~  I(ln_suitability*ln_price) + factor(year), index = "codigo", model = "within", effect = "individual", data = dta_r)
se3 <- coeftest(m3, vcov = vcovHC(m3, type = "HC2", cluster = "group")); se3

# Model 4: Restricted Sample, w/ controls
m4 <- plm(ln_conflict ~  I(ln_suitability*ln_price) + distance_hubs + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = "codigo", model = "within", effect = "individual", data = dta_r)
se4 <- coeftest(m4, vcov = vcovHC(m4, type = "HC2", cluster = "group")); se4

# Model 5: All Rural Unrest, w/o controls
m5 <- plm(ln_all ~ I(ln_suitability*ln_price) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m5) 
se5 <- coeftest(m5, vcov = vcovHC(m5, type = "HC2", cluster = "group")); se5

# Model 6: All Rural Unrest, w/ controls
m6 <- plm(ln_all ~ I(ln_suitability*ln_price) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m6) 
se6 <- coeftest(m6, vcov = vcovHC(m6, type = "HC2", cluster = "group")); se6

# Create table
stargazer(
  m1, m2, m3, m4, m5, m6,
  se = list(se1[,2], se2[,2], se3[,2], se4[,2], se5[,2], se6[,2]), 
  p = list(se1[,4], se2[,4], se3[,4], se4[,4], se5[,2], se6[,2]),
  omit = c("factor", "Constant", "distance_hubs", "untitled_share", "ln_pop"),
  omit.stat = c("adj.rsq", "ser", "f", "ll", "theta", "wald", "aic"),
  covariate.labels = c("Prices $\\times$ Suitability"), 
  add.lines = list(c("","","","","",""), 
                   c("Controls", rep(c("No", "Yes"), 3))),
  #c("Municipality FE", rep("Yes", 4)),
  #c("Year FE", rep("Yes", 4))),
  dep.var.labels.include = FALSE,
  model.names = FALSE,
  column.labels = c("\\text{\\normalfont Binary}","\\text{\\normalfont Binary}", 
                    "\\text{\\normalfont Restricted}", "\\text{\\normalfont Restricted}",
                    "\\text{\\normalfont All}", "\\text{\\normalfont All}"),
  df = FALSE, 
  multicolumn = TRUE, 
  type = "latex",
  star.cutoffs = c(0.1, 0.05, 0.01),
  star.char = c("*", "**", "***"),
  #float = TRUE, 
  #float.env = "sidewaystable",
  column.sep.width = ".5pt",
  font.size = "footnotesize",
  style = "ajps", 
  title = "Assessing Reporting Bias",
  label = "alt_meas",
  notes.append = FALSE,
  notes = c("\\parbox[t]{13.5}{\\textit{Note}: All models include municipality and year fixed effects. Standard errors clustered by municipality in parentheses. Models 1-2 are probit regression models. Models 3-6 are OLS models.}",
            "$^{*}p<.1$, $^{**}p<.05$, $^{***}p<.01$"),
  out = paste0(main,"1_writing/2023/tables/table_binposall.tex")
)

remove(m1, m2, m3, m4, m5, m6, se1, se2, se3, se4, se5, se6)


#### Table B.2: Lagged Agricultural Prices ----

# Model 1: Lagged Prices x Suitability, w/o controls
m1 <- plm(ln_conflict ~ I(ln_suitability*ln_price_l1) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m1) 
se1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC2", cluster = "group")); se1

# Model 2: Lagged Prices x Suitability, w/ controls
m2 <- plm(ln_conflict ~ I(ln_suitability*ln_price_l1) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m2) 
se2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC2", cluster = "group")); se2

# Model 3: Lagged Prices x Suitability x Settlements, w/o controls
m3 <- plm(ln_conflict ~ I(ln_suitability*ln_price_l1) + I(ln_price_l1*subsistence_farms_median) + I(ln_suitability*ln_price_l1*subsistence_farms_median) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m3) 
se3 <- coeftest(m3, vcov = vcovHC(m3, type = "HC2", cluster = "group")); se3

# Model 4: Lagged Prices x Suitability x Settlements, w/ controls
m4 <- plm(ln_conflict ~ I(ln_suitability*ln_price_l1) + I(ln_price_l1*subsistence_farms_median) + I(ln_suitability*ln_price_l1*subsistence_farms_median) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m4) 
se4 <- coeftest(m4, vcov = vcovHC(m4, type = "HC2", cluster = "group")); se4

# Model 5: Lagged Prices x Suitability x Committees, w/o controls
m5 <- plm(ln_conflict ~ I(ln_suitability*ln_price_l1) + I(ln_price_l1*orgs_locales_median) + I(ln_suitability*ln_price_l1*orgs_locales_median) + factor(year), index = c("codigo"), model = "within", effect = "individual", na.action = na.omit, data = dta); summary(m5) 
se5 <- coeftest(m5, vcov = vcovHC(m5, type = "HC2", cluster = "group")); se5

# Model 6: Lagged Prices x Suitability x Committees, w/ controls
m6 <- plm(ln_conflict ~ I(ln_suitability*ln_price_l1) + I(ln_price_l1*orgs_locales_median) + I(ln_suitability*ln_price_l1*orgs_locales_median) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", na.action = na.omit, data = dta); summary(m6) 
se6 <- coeftest(m6, vcov = vcovHC(m6, type = "HC2", cluster = "group")); se6

# Create table
stargazer(
  m1, m2, m3, m4, m5, m6, 
  se = list(se1[,2], se2[,2], se3[,2], se4[,2], se5[,2], se6[,2]), 
  p = list(se1[,4], se2[,4], se3[,4], se4[,4], se5[,4], se6[,4]),
  omit = c("factor", "Constant", "distance_hubs", "untitled_share", "ln_pop"),
  omit.stat = c("adj.rsq", "ser", "f"),
  covariate.labels = c("Prices$_{t-1}$ $\\times$ Suitability", 
                       "Prices$_{t-1}$ $\\times$ Settlements", 
                       "Prices$_{t-1}$ $\\times$ Suitability $\\times$ Settlements", 
                       "Prices$_{t-1}$ $\\times$ Committees", 
                       "Prices$_{t-1}$ $\\times$ Suitability $\\times$ Committees"),
  add.lines = list(c("","","","","",""), 
                   c("Controls", "No", "Yes", "No", "Yes", "No", "Yes")), 
  dep.var.labels.include = FALSE,
  model.names = FALSE,
  df = FALSE, 
  multicolumn = TRUE, 
  type = "latex",
  star.cutoffs = c(0.1, 0.05, 0.01),
  star.char = c("*", "**", "***"),
  #float = TRUE, 
  #float.env = "sidewaystable",
  column.sep.width = ".5pt",
  font.size = "scriptsize",
  style = "ajps", 
  label = "lagged",
  notes.append = FALSE,
  notes = c("\\parbox[t]{20cm}{\\textit{Note}: All models include municipality and year fixed effects. Standard errors clustered by municipality in parentheses. Prices are lagged by one year.}",
            "$^{*}p<.1$, $^{**}p<.05$, $^{***}p<.01$"),
  out = paste0(main,"1_writing/2023/tables/table_lagged.tex")
)


#### Figure B.1: Marginal Effect of Lagged Agricultural Prices on Peasant Resistance by Land Suitability ----

# Based on model 1
m1 <- plm(ln_conflict ~ ln_suitability*ln_price_l1 + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta)

# Covariance matrix for cluster standard errors
covMat1 <- vcovHC(m1, type = "HC2", cluster = "group")

# Plot
pdf(file = paste0(main, "1_writing/2023/figures/inter_lagged.pdf"), width = 5, height = 5)
interaction_plot_continuous(m1, effect="ln_price_l1", moderator="ln_suitability", interaction="ln_suitability:ln_price_l1", varcov=covMat1, minimum=0.75, maximum=1.4, incr="default", num_points = 50, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=T, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Lagged Prices")
dev.off()


#### Figure B.2: Marginal Effect of Lagged Prices on Peasant Resistance by Land Suitability, Subsistence Settlements, and Peasant Committees ----

# Based on model 3 & 5 (Table 1)
m3 <- plm(ln_conflict ~ ln_suitability*ln_price_l1 + ln_price_l1*subsistence_farms_median + ln_suitability*ln_price_l1*subsistence_farms_median + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta)
m5 <- plm(ln_conflict ~ ln_suitability*ln_price_l1 + ln_price_l1*orgs_locales_median + ln_suitability*ln_price_l1*orgs_locales_median + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta)

# Covariance matrices for cluster standard errors
covMat3 <- vcovHC(m3, type = "HC2", cluster = "group")
covMat5 <- vcovHC(m5, type = "HC2", cluster = "group")

# Plot (a): subsistence settlements
pdf(file = paste0(main, "1_writing/2023/figures/inter_subsist_lagged.pdf"), width = 5, height = 5)
interaction_plot_threeway_ylim(m3, effect="ln_price_l1", moderator1="ln_suitability", moderator2="subsistence_farms_median", varcov=covMat3, minimum=0.75, maximum=1.4, incr="default", num_points = 100, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=F, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Lagged Prices")
dev.off()
# Plot (b): peasant committees
pdf(file = paste0(main, "1_writing/2023/figures/inter_organiz_lagged.pdf"), width = 5, height = 5)
interaction_plot_threeway_ylim(m5, effect="ln_price_l1", moderator1="ln_suitability", moderator2="orgs_locales_median", varcov=covMat5, minimum=0.75, maximum=1.4, incr="default", num_points = 50, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=F, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Lagged Prices")
dev.off()


#### Table B.3: Alternative Splits of Organizational Resources ----

# Model 1: Prices x Suitability x Settlements (1Q), w/o controls
m1 <- plm(ln_conflict ~ I(ln_suitability*ln_price) + I(ln_price*subsistence_farms_1Q) + I(ln_suitability*ln_price*subsistence_farms_1Q) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m1) 
se1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC2", cluster = "group")); se1

# Model 2: Prices x Suitability x Settlements (1Q), w/ controls
m2 <- plm(ln_conflict ~ I(ln_suitability*ln_price) + I(ln_price*subsistence_farms_1Q) + I(ln_suitability*ln_price*subsistence_farms_1Q) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m2) 
se2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC2", cluster = "group")); se2

# Model 3: Prices x Suitability x Committees (1Q), w/o controls
m3 <- plm(ln_conflict ~ I(ln_suitability*ln_price) + I(ln_price*orgs_locales_1Q) + I(ln_suitability*ln_price*orgs_locales_1Q) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m3) 
se3 <- coeftest(m3, vcov = vcovHC(m3, type = "HC2", cluster = "group")); se3

# Model 4: Prices x Suitability x Committees (1Q), w/ controls
m4 <- plm(ln_conflict ~ I(ln_suitability*ln_price) + I(ln_price*orgs_locales_1Q) + I(ln_suitability*ln_price*orgs_locales_1Q) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m4) 
se4 <- coeftest(m4, vcov = vcovHC(m4, type = "HC2", cluster = "group")); se4

# Model 5: Prices x Suitability x Settlements (3Q), w/o controls
m5 <- plm(ln_conflict ~ I(ln_suitability*ln_price) + I(ln_price*subsistence_farms_3Q) + I(ln_suitability*ln_price*subsistence_farms_3Q) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m5) 
se5 <- coeftest(m5, vcov = vcovHC(m5, type = "HC2", cluster = "group")); se5

# Model 6: Prices x Suitability x Settlements (3Q), w/ controls
m6 <- plm(ln_conflict ~ I(ln_suitability*ln_price) + I(ln_price*subsistence_farms_3Q) + I(ln_suitability*ln_price*subsistence_farms_3Q) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m6) 
se6 <- coeftest(m6, vcov = vcovHC(m6, type = "HC2", cluster = "group")); se6

# Model 7: Prices x Suitability x Committees (3Q), w/o controls
m7 <- plm(ln_conflict ~ I(ln_suitability*ln_price) + I(ln_price*orgs_locales_3Q) + I(ln_suitability*ln_price*orgs_locales_3Q) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m7) 
se7 <- coeftest(m7, vcov = vcovHC(m7, type = "HC2", cluster = "group")); se7

# Model 8: Prices x Suitability x Committees (3Q), w/ controls
m8 <- plm(ln_conflict ~ I(ln_suitability*ln_price) + I(ln_price*orgs_locales_3Q) + I(ln_suitability*ln_price*orgs_locales_3Q) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m8) 
se8 <- coeftest(m8, vcov = vcovHC(m8, type = "HC2", cluster = "group")); se8

# Create table 
stargazer(
  m1, m2, m3, m4, m5, m6, m7, m8,
  se = list(se1[,2], se2[,2], se3[,2], se4[,2], se5[,2], se6[,2], se7[,2], se8[,2]), 
  p = list(se1[,4], se2[,4], se3[,4], se4[,4], se5[,4], se6[,4], se7[,4], se8[,4]),  
  covariate.labels = c("Prices $\\times$ Suitability",
                       "Prices $\\times$ Settlements (1Q)", 
                       "Prices $\\times$ Suitability $\\times$ Settlements (1Q)",
                       "Prices $\\times$ Committees (1Q)", 
                       "Prices $\\times$ Suitability $\\times$ Committees (1Q)",
                       "Prices $\\times$ Settlements (3Q)", 
                       "Prices $\\times$ Suitability $\\times$ Settlements (3Q)",
                       "Prices $\\times$ Committees (3Q)", 
                       "Prices $\\times$ Suitability $\\times$ Committees (3Q)" 
  ),
  add.lines = list(c("","","","","",""), 
                   c("Controls", rep(c("No", "Yes"), 4))),
  #c("Municipality FE", rep("Yes", 8)),
  #c("Year FE", rep("Yes", 8))),
  omit = c("factor", "Constant", "distance_hubs", "untitled_farmland", "ln_pop"),
  omit.stat = c("adj.rsq", "f"),
  dep.var.labels.include = FALSE,
  df = FALSE, 
  multicolumn = TRUE, 
  star.cutoffs = c(0.1, 0.05, 0.01),
  star.char = c("*", "**", "***"),
  #float = TRUE, 
  #float.env = "sidewaystable",
  #column.sep.width = "-1.5pt",
  font.size = "footnotesize",
  style = "ajps", 
  title = "Alternative Splits of Organizational Resources (1st and 3rd Quantiles)", 
  label = "alt_cutoffs",
  type = "latex",
  notes.append = FALSE,
  notes = c("\\parbox[t]{25.33cm}{\\textit{Note}: All models include municipality and year fixed effects. Standard errors clustered by municipality in parentheses. Models 1-4 are estimated on municipalities whose number of subsistence settlements and peasant committees are equal or less than the 1st quantile. Models 4-8 are estimated on municipalities whose number of subsistence settlements and peasant committes are greater than the 3rd quantile. The 1st and 3rd quantiles of subsistence settlements and peasant committees 2 and 8 and 0 and 5, respectively.}", 
            "$^{*}p<.1$, $^{**}p<.05$, $^{***}p<.01$"),
  out = paste0(main,"1_writing/2023/tables/table_cutoffs.tex")
)


#### Figure B.3: Marginal Effect of Agricultural Prices on Peasant Resistance by Land Suitability, Subsistence Settlements, and Peasant Committees (1st Quantile) ----

# Based on models 1 and 3
m1 <- plm(ln_conflict ~ ln_suitability*ln_price + ln_price*subsistence_farms_1Q + ln_suitability*ln_price*subsistence_farms_1Q + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta)
m3 <- plm(ln_conflict ~ ln_suitability*ln_price + ln_price*orgs_locales_1Q + ln_suitability*ln_price*orgs_locales_1Q + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta) 

# Covariance matrix for cluster standard errors
covMat1 <- vcovHC(m1, type = "HC2", cluster = "group")
covMat3 <- vcovHC(m3, type = "HC2", cluster = "group")

# Plot (a): subsistence settlements
pdf(file = paste0(main, "1_writing/2023/figures/inter_subsist_1Q.pdf"), width = 5, height = 5)
interaction_plot_threeway_ylim(m1, effect="ln_price", moderator1="ln_suitability", moderator2="subsistence_farms_1Q", varcov=covMat1, minimum=0.75, maximum=1.4, incr="default", num_points = 100, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=F, title=NULL, xlabel="Land suitability", ylabel="Marginal Effect of Agricultural Prices")
dev.off()
# Plot (b): peasant committees
pdf(file = paste0(main, "1_writing/2023/figures/inter_organiz_1Q.pdf"), width = 5, height = 5)
interaction_plot_threeway_ylim(m3, effect="ln_price", moderator1="ln_suitability", moderator2="orgs_locales_1Q", varcov=covMat3, minimum=0.75, maximum=1.4, incr="default", num_points = 50, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=F, title=NULL, xlabel="Land suitability", ylabel="Marginal Effect of Agricultural Prices")
dev.off()


#### Figure B.4: Marginal Effect of Agricultural Prices on Peasant Resistance by Land Suitability, Subsistence Settlements, and Peasant Committees (3rd Quantile) ----

# Based on models 5 and 7 
m5 <- plm(ln_conflict ~ ln_suitability*ln_price + ln_price*subsistence_farms_3Q + ln_suitability*ln_price*subsistence_farms_3Q + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta)
m7 <- plm(ln_conflict ~ ln_suitability*ln_price + ln_price*orgs_locales_3Q + ln_suitability*ln_price*orgs_locales_3Q + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta) 

# Covariance matrix for cluster standard errors
covMat5 <- vcovHC(m5, type = "HC2", cluster = "group")
covMat7 <- vcovHC(m7, type = "HC2", cluster = "group")

# Plot (a): subsistence settlements
pdf(file = paste0(main, "1_writing/2023/figures/inter_subsist_3Q.pdf"), width = 5, height = 5)
interaction_plot_threeway_ylim(m5, effect="ln_price", moderator1="ln_suitability", moderator2="subsistence_farms_3Q", varcov=covMat5, minimum=0.75, maximum=1.4, incr="default", num_points = 100, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=F, title=NULL, xlabel="Land suitability", ylabel="Marginal Effect of Agricultural Prices")
dev.off()
# Plot (b): peasant committees
pdf(file = paste0(main, "1_writing/2023/figures/inter_organiz_3Q.pdf"), width = 5, height = 5)
interaction_plot_threeway_ylim(m7, effect="ln_price", moderator1="ln_suitability", moderator2="orgs_locales_3Q", varcov=covMat7, minimum=0.75, maximum=1.4, incr="default", num_points = 50, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=F, title=NULL, xlabel="Land suitability", ylabel="Marginal Effect of Agricultural Prices")
dev.off()

remove(m1, m2, m3, m4, m5, m6, m7, m8, se1, se2, se3, se4, se5, se6, se7, se8, covMat1, covMat3, covMat5, covMat7, dta_r, zeros)


#### Table B.4: Alternative Base Regression Models ----

# Model 1: Department FEs x Year FEs, w/o controls
m1 <- lm(ln_conflict ~ priceXsuit + factor(dept_year) + factor(year), data = dta); summary(m1) 
se1 <- coeftest(m1, vcov = vcovCL(m1, cluster = ~codigo)); se1

# Model 2: Department FEs x Year FEs, w/ controls
m2 <- lm(ln_conflict ~ priceXsuit + factor(dept_year) + factor(year)*distance_hubs + factor(year)*untitled_share, data = dta); summary(m2) 
se2 <- coeftest(m2, vcov = vcovCL(m2, cluster = ~codigo)); se2

# Model 3: Random Effects, w/o controls
m3 <- plm(ln_conflict ~ ln_price + ln_suitability + priceXsuit, index = c("codigo", "year"), model = "random", effect = "twoway", data = dta); summary(m3)
se3 <- coeftest(m3, vcov = vcovHC(m3, type = "HC2", cluster = "group")); se3

# Model 4: Random Effects, w/ controls
m4 <- plm(ln_conflict ~ ln_price + ln_suitability + priceXsuit + distance_hubs + untitled_share + ln_pop, index = c("codigo", "year"), model = "random", effect = "twoway", data = dta); summary(m4)
se4 <- coeftest(m4, vcov = vcovHC(m4, type = "HC2", cluster = "group")); se4

# Model 5: Lagged Dependent Variable, w/o controls
m5 <- lm(ln_conflict ~ ln_price + ln_suitability + priceXsuit + ln_conflict_t1 + ln_conflict_t2 + ln_conflict_t3, data = dta); summary(m5)
se5 <- pcse(m5, groupN = dta$codigo, groupT = dta$year); se5

# Model 6: Lagged Dependent Variable, w/ controls
m6 <- lm(ln_conflict ~ ln_price + ln_suitability + priceXsuit + ln_conflict_t1 + ln_conflict_t2 + ln_conflict_t3 + distance_hubs + untitled_share + ln_pop, subset = !(is.na(untitled_share)), data = dta); summary(m6)
se6 <- pcse(m6, groupN = subset(dta, !(is.na(untitled_share)))$codigo, groupT = subset(dta, !(is.na(untitled_share)))$year); se6

# Model 7: Tobit, w/o controls
m7 <- tobit(ln_conflict ~ ln_price + ln_suitability + priceXsuit, left = 0, robust = TRUE, data = dta); summary(m7)

# Model 8: Tobit, w/ controls
m8 <- tobit(ln_conflict ~ ln_price + ln_suitability + priceXsuit + distance_hubs + untitled_share + ln_pop, left = 0, robust = TRUE, data = dta); summary(m8)

# Model 9: Negative Binomial, w/o controls
m9 <- glm.nb(n_conflict ~ priceXsuit + factor(codigo) + factor(year), data = dta); summary(m9)
se9 <- coeftest(m9, vcov = vcovHC(m9, type = "HC2", cluster = "group")); se9

# Model 10: Negative Binomial, w/ controls
m10 <- glm(n_conflict ~ priceXsuit + factor(year)*distance_hubs + factor(year)*untitled_share + ln_pop + factor(codigo), data = dta); summary(m10)
se10 <- coeftest(m10, vcov = vcovHC(m10, type = "HC2", cluster = "group")); se10

# Create table 
stargazer(
  m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, 
  se = list(se1[,2], se2[,2], se3[,2], se4[,2], se5[["pcse"]], se6[["pcse"]], NULL, NULL, se9[,2], se10[,2]), 
  p = list(se1[,4], se2[,4], se3[,4], se4[,4], se5[["pval"]], se6[["pval"]], NULL, NULL, se9[,4], se10[,4]),  
  omit = c("factor", "Constant", "distance_hubs", "untitled_share", "ln_pop",
           "ln_conflict_t1", "ln_conflict_t2", "ln_conflict_t3"),
  covariate.labels = c("Prices", 
                       "Suitability", 
                       "Prices $\\times$ Suitability"),
  add.lines = list(c("","","","","","","","","",""),
                   c("Controls", rep(c("No", "Yes"), 5)),
                   c("Department FE $\\times$ Year FE", rep("Yes", 2), rep("No", 8)),
                   c("Municipality FE", rep("No", 8), rep("Yes", 2)),
                   #c("Department FE", rep("Yes", 2), rep("No", 8)),
                   c("Year FE", rep("Yes", 2), rep("No", 6), rep("Yes", 2)), 
                   c("Lagged DV", rep("No", 4), rep("Yes", 2), rep("No", 4))),
  omit.stat = c("adj.rsq", "ser", "f", "ll", "theta", "wald", "aic"),
  dep.var.labels.include = FALSE,
  model.names = FALSE,
  column.labels = c("$FE$","$FE$","$RE$","$RE$","$AR$","$AR$","$Tobit$","$Tobit$","$NB$","$NB$"),
  df = FALSE, 
  multicolumn = TRUE, 
  star.cutoffs = c(0.1, 0.05, 0.01),
  star.char = c("*", "**", "***"),
  float = TRUE, 
  float.env = "sidewaystable",
  font.size = "footnotesize",
  style = "ajps", 
  title = "Alternative Base Regression Models", 
  label = "alt_models",
  type = "latex",
  notes.append = FALSE,
  notes = c("\\parbox[t]{23.5cm}{\\textit{Note}: Standard errors clustered by municipality (models 1-4 and 9-10), panel corrected (models 5-6), and Huber-White robust (models 7-8) in parentheses. Constants estimated but not reported. The unit of analysis is the municipality-year. Models 1-2 include interactive fixed effects between department and year, and year fixed effects. Models 3-4 include two-way random effects. Models 5-6 include auto-regressive models with three lagged terms of the dependent variable. Models 7-8 are tobit regression models for censored data. Models 9-10 are hybrid, fixed-effects negative binomial models for event-count data. Unit-level fixed effects are omitted in models 5-6 because OLS estimates are biased in models with a lagged dependent variable and fixed effects \\citep{Nickell1981}. They are also omitted in models 7-8 because of the incidental parameter problem, which yields biased and inconsistent estimates in tobit models \\citep{Greene2004}. Hybrid negative binomial estimators are preferred because conditional maximum likelihood approaches for negative binomial are not true fixed-effects models \\citep{AllisonWaterman2002}.}",
            "$^{*}p<.1$, $^{**}p<.05$, $^{***}p<.01$"),
  out = paste0(main,"1_writing/2023/tables/table_alt_models.tex")
)

remove(m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, se1, se2, se3, se4, se5, se6, se9, se10)


#### Figure B.5: Binning Estimates of the Marginal Effect of Agricultural Prices on Peasant Resistance by Land Suitability ----

# Binning two-way
binning <- interflex(
  "binning", dta, Y = "ln_conflict", D = "ln_price", X = "ln_suitability", FE = NULL, cl = "codigo", na.rm = TRUE, cex.lab = 1, 
  xlab = "Land Suitability", ylab = "Marginal Effect of Agricultural Prices",
  theme.bw = TRUE, show.grid = FALSE)

# Plot
binning.plot <- binning$figure
pdf(file = paste0(main, "1_writing/2023/figures/inter_binning_twoway.pdf"), width = 6, height = 6)
binning.plot
dev.off()


#### Table B.5: Binning Estimates by Land Suitability ----

# Estimates of binning plots
est.bin <- as.data.frame(binning[["est.bin"]])
rownames(est.bin) <- substr(rownames(est.bin), 16, nchar(rownames(est.bin)))
colnames(est.bin)[1] <- "x0"
colnames(est.bin)[2] <- "coef"
colnames(est.bin)[3] <- "sd"
colnames(est.bin)[4] <- "CI_lower"
colnames(est.bin)[5] <- "CI_upper"

# Table
binning.est <- xtable(rbind(
  c(paste0("L: ", rownames(est.bin)[1]), 
    paste0("M: ", rownames(est.bin)[2]),
    paste0("H: ", rownames(est.bin)[3])),
  round(est.bin$coef, 3), 
  c(paste0("[", round(est.bin$CI_lower[1], 3), ",", round(est.bin$CI_upper[1], 3), "]"), 
    paste0("[", round(est.bin$CI_lower[2], 3), ",", round(est.bin$CI_upper[2], 3), "]"), 
    paste0("[", round(est.bin$CI_lower[3], 3), ",", round(est.bin$CI_upper[3], 3), "]")), 
  #c(round(binning[["model.binning"]][["rpval"]][["D.G.1"]], 3), 
  #  round(binning[["model.binning"]][["rpval"]][["D.G.2"]], 3),
  #  round(binning[["model.binning"]][["rpval"]][["D.G.3"]], 3)),
  c("$\\alpha_{2} = \\alpha_{1}$", "$\\alpha_{3} = \\alpha_{2}$", "$\\alpha_{1} = \\alpha_{3}$")))

row.names(binning.est) <- c('Bins', 'Estimate', '95\\% CI', #'$p$-value', 
                            "$H_{0}$")
colnames(binning.est) <- c("", "\\textit{Land Suitability}", "")
binning.est

addtorow <- list()
addtorow$pos <- list(4)
addtorow$command <- c("\\hline \\addlinespace[0.5em] \\multicolumn{4}{l}{\\parbox[t]{12cm}{Marginal effect of agricultural prices on peasant resistance. Based on Figure \\ref{binning_twoway}. L, M, and H refer to the first, second, and third tercile of land suitability.}} \\\\")

sink(paste0(main,"1_writing/2023/tables/table_binning_2014.tex"))
print(xtable(binning.est, type = "latex"),
      booktabs = TRUE, 
      include.colnames = TRUE,
      sanitize.text.function = function(x){x},
      hline.after = c(-1, 0),
      size = "normalsize", 
      floating = FALSE, 
      add.to.row = addtorow)
sink()


#### Figure B.6: Binning Estimates of the Marginal Effect of Agricultural Prices on Peasant Resistance by Land Suitability, Subsistence Settlements, and Peasant Committees ----

## Plot (a): subsistence settlements
binning.subsist0 <- interflex(
  "binning", subset(dta, subsistence_farms_median == 0),
  Y = "ln_conflict", D = "ln_price", X = "ln_suitability", Z = NULL, FE = NULL, cl = "codigo", na.rm = TRUE, 
  xlab = "Land Suitability", ylab = "Marginal Effect of Agricultural Prices",
  theme.bw = TRUE, show.grid = FALSE)

binning.subsist0.plot <- binning.subsist0$figure
pdf(file = paste0(main, "1_writing/2023/figures/inter_binning_subsist0.pdf"), width = 6, height = 6)
binning.subsist0.plot
dev.off()

binning.subsist1 <- interflex(
  "binning", subset(dta, subsistence_farms_median == 1),
  Y = "ln_conflict", D = "ln_price", X = "ln_suitability", Z = NULL, FE = NULL, cl = "codigo", na.rm = TRUE, 
  xlab = "Land Suitability", ylab = "Marginal Effect of Agricultural Prices",
  theme.bw = TRUE, show.grid = FALSE)

binning.subsist1.plot <- binning.subsist1$figure
pdf(file = paste0(main, "1_writing/2023/figures/inter_binning_subsist1.pdf"), width = 6, height = 6)
binning.subsist1.plot
dev.off()

## Plot (b): peasant committees
binning.organiz0 <- interflex(
  "binning", subset(dta, orgs_locales_median == 0),
  Y = "ln_conflict", D = "ln_price", X = "ln_suitability", Z = "year.factor", FE = NULL, cl = "codigo", na.rm = TRUE, 
  xlab = "Land Suitability", ylab = "Marginal Effect of Agricultural Prices",
  theme.bw = TRUE, show.grid = FALSE)

binning.organiz0.plot <- binning.organiz0$figure + 
  pdf(file = paste0(main, "1_writing/2023/figures/inter_binning_organiz0.pdf"), width = 6, height = 6)
binning.organiz0.plot
dev.off()

binning.organiz1 <- interflex(
  "binning", subset(dta, orgs_locales_median == 1),
  Y = "ln_conflict", D = "ln_price", X = "ln_suitability", Z = "year.factor", FE = NULL, cl = "codigo", na.rm = TRUE, 
  xlab = "Land Suitability", ylab = "Marginal Effect of Agricultural Prices",
  theme.bw = TRUE, show.grid = FALSE)

binning.organiz1.plot <- binning.organiz1$figure
pdf(file = paste0(main, "1_writing/2023/figures/inter_binning_organiz1.pdf"), width = 6, height = 6)
binning.organiz1.plot
dev.off()


#### Table B.6: Binning Estimates by Land Suitability and Subsistence Agriculture ----

# Estimates of binning plots
est.bin.sub0 <- as.data.frame(binning.subsist0[["est.bin"]])
rownames(est.bin.sub0) <- substr(rownames(est.bin.sub0), 16, nchar(rownames(est.bin.sub0)))
colnames(est.bin.sub0)[1] <- "x0"
colnames(est.bin.sub0)[2] <- "coef"
colnames(est.bin.sub0)[3] <- "sd"
colnames(est.bin.sub0)[4] <- "CI_lower"
colnames(est.bin.sub0)[5] <- "CI_upper"

est.bin.sub1 <- as.data.frame(binning.subsist1[["est.bin"]])
rownames(est.bin.sub1) <- substr(rownames(est.bin.sub1), 16, nchar(rownames(est.bin.sub1)))
colnames(est.bin.sub1)[1] <- "x0"
colnames(est.bin.sub1)[2] <- "coef"
colnames(est.bin.sub1)[3] <- "sd"
colnames(est.bin.sub1)[4] <- "CI_lower"
colnames(est.bin.sub1)[5] <- "CI_upper"

# Table
binning.est.subsist <- xtable(rbind(
  c("", "", ""),
  # W/o settlements:
  c("$\\leq$ \\textbf{\\textit{med}(Settlements) \\medskip}", "", ""),
  c(paste0("L: ", rownames(est.bin.sub0)[1]), 
    paste0("M: ", rownames(est.bin.sub0)[2]),
    paste0("H: ", rownames(est.bin.sub0)[3])),
  round(est.bin.sub0$coef, 3), 
  c(paste0("[", round(est.bin.sub0$CI_lower[1], 3), ",", round(est.bin.sub0$CI_upper[1], 3), "]"), 
    paste0("[", round(est.bin.sub0$CI_lower[2], 3), ",", round(est.bin.sub0$CI_upper[2], 3), "]"), 
    paste0("[", round(est.bin.sub0$CI_lower[3], 3), ",", round(est.bin.sub0$CI_upper[3], 3), "]")), 
  #c(round(binning.subsist0[["model.binning"]][["rpval"]][["D.G.1"]], 3), 
  #  round(binning.subsist0[["model.binning"]][["rpval"]][["D.G.2"]], 3),
  #  round(binning.subsist0[["model.binning"]][["rpval"]][["D.G.3"]], 3)),
  c("$\\alpha_{2} = \\alpha_{1}$", "$\\alpha_{3} = \\alpha_{2}$", "$\\alpha_{1} = \\alpha_{3}$"),
  c("", "", ""),
  # W/ settlements:
  c("$>$ \\textbf{\\textit{med}(Settlements) \\medskip}", "", ""),
  c(paste0("L: ", rownames(est.bin.sub1)[1]), 
    paste0("M: ", rownames(est.bin.sub1)[2]),
    paste0("H: ", rownames(est.bin.sub1)[3])),
  round(est.bin.sub1$coef, 3), 
  c(paste0("[", round(est.bin.sub1$CI_lower[1], 3), ",", round(est.bin.sub1$CI_upper[1], 3), "]"), 
    paste0("[", round(est.bin.sub1$CI_lower[2], 3), ",", round(est.bin.sub1$CI_upper[2], 3), "]"), 
    paste0("[", round(est.bin.sub1$CI_lower[3], 3), ",", round(est.bin.sub1$CI_upper[3], 3), "]")), 
  #c(round(binning.subsist1[["model.binning"]][["rpval"]][["D.G.1"]], 3), 
  #  round(binning.subsist1[["model.binning"]][["rpval"]][["D.G.2"]], 3),
  #  round(binning.subsist1[["model.binning"]][["rpval"]][["D.G.3"]], 3)),
  c("$\\alpha_{2} = \\alpha_{1}$", "$\\alpha_{3} = \\alpha_{2}$", "$\\alpha_{1} = \\alpha_{3}$")
))

binning.est.subsist <- as.matrix(binning.est.subsist)
rownames(binning.est.subsist) <- c("", "\\textbf{(a)}", 'Bins', 'Estimate', '95\\% CI', "$H_{0}$", "  ", "\\textbf{(b)}",  'Bins ', 'Estimate ', '95\\% CI ', "$H_{0}$ ")
# '\\textit{p-}value', '\\textit{p-}value '
colnames(binning.est.subsist) <- c(" ", "\\textit{Land Suitability}", "  ")
binning.est.subsist

addtorow <- list()
addtorow$pos <- list(12)
addtorow$command <- c("\\hline \\addlinespace[0.5em] \\multicolumn{4}{l}{\\parbox[t]{12cm}{Marginal effect of agricultural prices on peasant resistance. Based on Figure \\ref{binning_threeway}. L, M, and H refer to the first, second, and third tercile of land suitability.}} \\\\")

sink(paste0(main,"1_writing/2023/tables/table_binning_subsist_2014.tex"))
print(xtable(binning.est.subsist, type = "latex"),
      booktabs = TRUE, 
      include.colnames = TRUE,
      sanitize.text.function = function(x){x},
      sanitize.rownames.function = identity,
      hline.after = c(-1, 0),
      size = "normalsize", 
      floating = FALSE, 
      add.to.row = addtorow
)
sink()


#### Table B.7: Binning Estimates by Land Suitability and Organizational Capacities ----

# Estimates of binning plots
est.bin.org0 <- as.data.frame(binning.organiz0[["est.bin"]])
rownames(est.bin.org0) <- substr(rownames(est.bin.org0), 16, nchar(rownames(est.bin.org0)))
colnames(est.bin.org0)[1] <- "x0"
colnames(est.bin.org0)[2] <- "coef"
colnames(est.bin.org0)[3] <- "sd"
colnames(est.bin.org0)[4] <- "CI_lower"
colnames(est.bin.org0)[5] <- "CI_upper"

est.bin.org1 <- as.data.frame(binning.organiz1[["est.bin"]])
rownames(est.bin.org1) <- substr(rownames(est.bin.org1), 16, nchar(rownames(est.bin.org1)))
colnames(est.bin.org1)[1] <- "x0"
colnames(est.bin.org1)[2] <- "coef"
colnames(est.bin.org1)[3] <- "sd"
colnames(est.bin.org1)[4] <- "CI_lower"
colnames(est.bin.org1)[5] <- "CI_upper"

# Table
binning.est.organiz <- xtable(rbind(
  c("", "", ""),
  # W/o committees:
  c("$\\leq$ \\textbf{\\textit{med}(Committees) \\medskip}", "", ""),
  c(paste0("L: ", rownames(est.bin.org0)[1]), 
    paste0("M: ", rownames(est.bin.org0)[2]),
    paste0("H: ", rownames(est.bin.org0)[3])),
  round(est.bin.org0$coef, 3), 
  c(paste0("[", round(est.bin.org0$CI_lower[1], 3), ",", round(est.bin.org0$CI_upper[1], 3), "]"), 
    paste0("[", round(est.bin.org0$CI_lower[2], 3), ",", round(est.bin.org0$CI_upper[2], 3), "]"), 
    paste0("[", round(est.bin.org0$CI_lower[3], 3), ",", round(est.bin.org0$CI_upper[3], 3), "]")), 
  #c(round(binning.organiz0[["model.binning"]][["rpval"]][["D.G.1"]], 3), 
  #  round(binning.organiz0[["model.binning"]][["rpval"]][["D.G.2"]], 3),
  #  round(binning.organiz0[["model.binning"]][["rpval"]][["D.G.3"]], 3)),
  c("$\\alpha_{2} = \\alpha_{1}$", "$\\alpha_{3} = \\alpha_{2}$", "$\\alpha_{1} = \\alpha_{3}$"),
  c("", "", ""),
  # W/ committees:
  c("$>$ \\textbf{\\textit{med}(Committees) \\medskip}", "", ""),
  c(paste0("L: ", rownames(est.bin.org1)[1]), 
    paste0("M: ", rownames(est.bin.org1)[2]),
    paste0("H: ", rownames(est.bin.org1)[3])),
  round(est.bin.org1$coef, 3), 
  c(paste0("[", round(est.bin.org1$CI_lower[1], 3), ",", round(est.bin.org1$CI_upper[1], 3), "]"), 
    paste0("[", round(est.bin.org1$CI_lower[2], 3), ",", round(est.bin.org1$CI_upper[2], 3), "]"), 
    paste0("[", round(est.bin.org1$CI_lower[3], 3), ",", round(est.bin.org1$CI_upper[3], 3), "]")), 
  #c(round(binning.organiz1[["model.binning"]][["rpval"]][["D.G.1"]], 3), 
  #  round(binning.organiz1[["model.binning"]][["rpval"]][["D.G.2"]], 3),
  #  round(binning.organiz1[["model.binning"]][["rpval"]][["D.G.3"]], 3)),
  c("$\\alpha_{2} = \\alpha_{1}$", "$\\alpha_{3} = \\alpha_{2}$", "$\\alpha_{1} = \\alpha_{3}$")
))

binning.est.organiz <- as.matrix(binning.est.organiz)
rownames(binning.est.organiz) <- c("", "\\textbf{(a)}", 'Bins', 'Estimate', '95\\% CI', "$H_{0}$", "  ", "\\textbf{(b)}",  'Bins ', 'Estimate ', '95\\% CI ', "$H_{0}$ ")
# '\\textit{p-}value', '\\textit{p-}value '
colnames(binning.est.organiz) <- c(" ", "\\textit{Land Suitability}", "  ")
binning.est.organiz

addtorow <- list()
addtorow$pos <- list(12)
addtorow$command <- c("\\hline \\addlinespace[0.5em] \\multicolumn{4}{l}{\\parbox[t]{12cm}{Marginal effect of agricultural prices on peasant resistance. Based on Figure \\ref{binning_threeway}, L, M, and H refer to the first, second, and third tercile of land suitability.}} \\\\")

sink(paste0(main,"1_writing/2023/tables/table_binning_organiz_2014.tex"))
print(xtable(binning.est.organiz, type = "latex"),
      booktabs = TRUE, 
      include.colnames = TRUE,
      sanitize.text.function = function(x){x},
      sanitize.rownames.function = identity,
      hline.after = c(-1, 0),
      size = "normalsize", 
      floating = FALSE, 
      add.to.row = addtorow
)
sink()

remove(binning, binning.plot, binning.est, addtorow, est.bin.sub0, est.bin.sub1, est.bin.org0, est.bin.org1,
       binning.subsist0, binning.subsist0.plot, binning.subsist1, binning.subsist1.plot, binning.est.subsist,
       binning.organiz0, binning.organiz0.plot, binning.organiz1, binning.organiz1.plot, binning.est.organiz)

remove(m1, m2, m3, m4, m5, m6, m7, m8, se1, se2, se3, se4, se5, se6, s7, se8)


